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1.  Introduction 


The  detection  and  identification  of  landmines  and  unexploded  ordnance  (UXO)  is  a  difficult 
problem  that  has  continuously  challenged  our  military’s  scientific  and  technological  expertise. 

Of  the  various  mitigation  efforts  the  U.S.  Army  has  undertaken  over  the  years  in  response  to 
these  threats,  sensing  with  radio  waves  has  garnered  significant  attention — and  also  has 
gradually  become  the  most  prominent  technique.  In  this  method,  a  transmitter  illuminates  a  scene 
with  electromagnetic  signals  and  then  the  reflected — or  scattered — signals  are  recorded  and 
processed  to  identify  the  presence  and  locations  of  potential  targets.  As  a  target’s  signature  is 
often  distorted  or  obscured  by  the  responses  from  the  environment,  the  ability  to  model  and 
simulate  the  target’s  electromagnetic  interactions  with  its  surroundings  is  therefore  critical  in  the 
development  of  high  fidelity  detection  systems. 

An  important  sub-problem  in  electromagnetic-oriented  sensing  of  both  surface  and  underground 
targets  is  the  evaluation  and  understanding  of  ground  surface  scattering  behavior.  Although 
many  phenomenological  studies  for  natural  surfaces  have  been  carried  out  within  the  radio 
frequency  (RF)  remote  sensing  community  on  the  inference  of  physical  parameters  (such  as  soil 
composition  and  moisture  content)  from  the  surfaces’  polarimetric  and  spectral  scattering 
signature,  the  emphasis  of  contemporary  investigations  has  been  primarily  on  applications 
related  to  airborne  or  spacebome  technologies.  As  such,  the  standard  radar  operational  modality 
of  interest  is  often  confined  to  non-grazing  observation  angles  while  the  operational  frequency  is 
at  L-band  and  above.  Recently,  there  has  been  considerable  interest  in  the  development  of 
ground-based  mobile  sensing  platforms  for  standoff  detection  and  identification  of  in-road  and 
roadside  threats  (1-4).  In  these  low-to-ground  systems,  as  the  radio  wave  propagation  paths 
defining  the  electromagnetic  interactions  between  the  radar  transceiver  and  targets  adhere  to  the 
ground,  existing  ground  surface  scattering  models  must  be  supplemented,  or  extended,  to  include 
low  grazing  angle  effects.  Accordingly,  a  full-wave  electromagnetic  simulation  approach  is 
proposed  in  this  work  to  estimate  rough  surface  background  clutter  as  relevant  to  performance 
prediction  for  the  forward-looking  imaging  radar  developed  at  the  U.S.  Army  Research 
Laboratory  (ARL). 

The  low-frequency,  ultra-wideband,  synchronous  impulse  reconstruction  (SIRE)  radar  testbed 
(7)  designed  by  ARL  is  a  vehicle-based  system  operating  at  the  nominal  frequency  band  of 
500-2500  MHz  with  the  forward-looking  coverage  angle  6inc  (figure  1)  spanning  approximately 
5°-15°  with  respect  to  the  horizon.  Although  the  performance  of  the  radar  can  be  limited  by 
either  the  -target-to-system  noise”  or  the  -target-to-background”  ratio,  the  subject  of  interest  in 
the  current  study  is  the  latter  quantity  as  it  directly  defines  the  theoretical  physics-based 
detectability  limits  of  the  radar.  A  salient  feature  associated  with  scattering  from  surface  targets 
at  grazing  angles  is  that  the  signal  illuminating  the  target,  as  well  as  the  signal  subsequently  re- 
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radiated  by  the  target,  both  experience  substantial  attenuations  owing  to  the  near  cancellation  of 
their  direct  and  ground-reflected  wavefronts.  As  a  result,  the  backscattered  signal  from  the 
ground  itself  may  turn  out  to  be  strong  enough  to  act  as  an  interference  component  in  masking 
the  target  response.  Explicitly,  it  is  the  relative  strength  of  these  two  signals  that  needs  to  be 
carefully  studied,  especially  for  a  small  discrete  target. 


Figure  1.  Electromagnetic  wave  illuminating  a  rough  surface  at  grazing  angle  9inc. 

The  strength  and  pattern  of  the  backscattered  signal  from  the  ground  are  complicated  functions 
of  the  radar  geometry  (i.e.,  signal  incidence  and  observation  angles),  transmission  and  reception 
polarizations,  frequency,  as  well  as  the  electrical  and  physical  properties  of  the  ground.  The 
unique  properties,  or  complications,  of  near-earth  propagation  such  as  surface  wave  propagation, 
non-plane  wave  propagation,  and  higher  order  reflection  and  diffraction  phenomena  pose 
additional  constraints  that  often  beset  the  validity  of  classical  analytical  ray-tracing  and  physical 
optics  techniques — and  their  heuristic  extensions.  Another  issue  in  need  of  further  investigation 
is  how  these  special  grazing  effects  can  be  exploited  for  target  signature  extraction  and  image 
clutter  suppression  pertaining  to  radar  applications.  Alas,  few  comprehensive  experimental  data 
sets  exist  for  characterizing  ground  clutter  at  the  frequency  band  stated  earlier.  In  open  literature, 
for  example,  area-wide  land  clutter  data  can  be  found  in  references  5  and  6  at  five  frequency 
bands  (very  high  frequency  [VHF],  ultra  high  frequency  [UHF],  and  L-,  S-,  and  A-bands),  but 
these  data  are  only  relevant  for  a  radar  with  a  field  of  view  spanning  over  many  kilometers  of 
terrain.  The  empirically  based  formula  presented  in  reference  7  (derived  from  measurements  in 
reference  8)  for  estimating  background  clutter  is  only  valid  at  the  higher  frequency  bands.  Fully 
polarimetric  measurements  of  soil  have  been  done  by  Oh  et  al.  (9)  at  L-,  C-,  and  X-bands  for  a 
limited  set  of  soil  parameters.  Overall,  there  is  a  scarcity  of  experimental  ground  scattering  data 
at  grazing  angles  and  low  frequencies. 
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Theoretical  studies  of  scattering  and  propagation  in  the  presence  of  an  undulating  interface 
separating  two  media  (e.g.,  air/ground)  have  been  undertaken  by  various  researchers.  Commonly 
employed  analytical  approaches  for  computing  electromagnetic  surface  scattering  include  the 
small  perturbation  method  (SPM)  (10)  and  the  Kirchhoff  approximation  (KA)  (11,  12).  The 
validity  of  these  theoretical  techniques  is  dependent  upon  the  physical  statistics  of  the  surface 
profile.  Specifically,  SPM  is  only  applicable  for  surfaces  with  root  mean  square  (rms)  height 
variations  much  smaller  than  the  wavelength,  whereas  KA  is  intended  for  surfaces  with 
undulations  that  are  large  compared  to  the  wavelength  (as  long  as  the  average  radius  of  curvature 
of  the  profile  is  also  large).  In  addition,  these  methods  do  not  produce  meaningful  results  when 
applied  to  modeling  scattering  and  propagation  phenomena  at  near-grazing  angles.  Over  the 
years,  much  effort  has  been  devoted  to  supplementing  perturbation  and  KA  scattering  models  to 
extend  their  range  of  validity.  Works  in  this  endeavor  include  high-order  SPM  (13,  14),  small 
slope  approximation  (15,  16),  phase  perturbation  method  (17),  unified  perturbation  expansion 
(18),  and  the  polarization  current-based  perturbation  expansion  (19,  20);  however,  the  emphasis 
of  these  studies  is  not  on  the  grazing  configuration.  Consequently,  closed- form  analytical 
expressions  for  the  radar  scattering  coefficient  derived  from  the  aforementioned  techniques  may 
be  of  limited  use  for  addressing  the  specific  needs  of  the  current  problem  of  interest. 

In  view  of  the  lack  of  empirical  data  and  the  deficiencies  of  analytical  treatments,  numerical 
simulations  of  surface  scattering  effects  based  on  direct  field  solvers — such  as  integral  equation 
and  finite  difference  methods — have  been  put  forth  in  previous  works  (21)  (and  references 
therein).  A  survey  of  these  studies  shows  that  existing  integral  equation  solvers  are  rather 
unattractive  for  the  three-dimensional  (3-D)  grazing  problem  and  very  few  of  the  proposed  finite 
difference  methods  focus  on  the  grazing  scenario  at  all.  The  twofold  objective  of  the  current 
study  is  discussed  as  follows.  First,  in  sections  2  and  3,  the  development  of  the  finite-difference 
time-domain  (FDTD)  technique  based  on  the  finite-extent  perturbation  surface  model  is 
outlined — including  a  detailed  investigation  of  the  validity  of  the  method  obtained  through 
comparisons  with  a  high-order  surface  integral  equation  (SIE)  solution  and  with  published 
measurement  data.  Then,  in  section  4,  a  demonstration  of  the  effects  of  ground  clutter  on  target 
imaging  generated  by  the  time-reversal  technique  is  presented  with  the  simulation  of  a  large 
terrain  scene  consisting  of  targets  in  a  rough  ground  environment.  In  sum,  this  study  illustrates 
the  practicality  of  the  chosen  full-wave  simulation  strategy  for  the  emulation  of  forward-looking 
radar  operation  and  imaging. 


2.  Validation  of  FDTD  Simulation  of  Rough  Surface  Scattering  in  2-D 


The  purpose  of  this  section  is  to  determine  the  accuracy  (and  limitations)  of  the  proposed  FDTD 
algorithm  for  various  grazing  incidence  angles  and  terrain  surface  parameters  by  using  a  high- 
order  SIE-derived  solution  as  the  reference  solution.  Motivations  are  also  given  to  illustrate  why 
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the  employed  approach  may  be  more  appropriate  than  the  SIE  solver  for  the  current  work  of 
interest.  To  this  end,  the  SIE  solver  is  discussed  first,  followed  by  the  FDTD.  The  two- 
dimensional  (2-D)  problem  is  considered  at  the  moment  to  reduce  the  computational  complexity 
of  the  validation  process. 

2.1  SIE  Solver 

The  high-order  accurate  integral  equation  solver  originally  developed  in  references  22  and  23  for 
long-distance,  near-ground  wave  propagation  applications  is  modified  to  characterize  scattering 
from  an  irregular  terrain  profile.  The  solver  discretizes  the  Poggio,  Miller,  Chang,  Harrington, 
Wu,  and  Tsai  (PMCHWT)  ( 24-26)  combined- field  surface  integral  formulation  using  the  regular 
locally  corrected  Nystrom  (LCN)  method  (27).  In  essence,  the  following  boundary  conditions  are 
observed  at  the  one-dimensional  (1-D)  ground  surface  in  relating  the  incident  electric  and 
magnetic  fields  to  the  induced  electric  and  magnetic  surface  currents: 

n  x  Emc  =  -h  x  [(jhL,  +  7j2L2  )J(t)  -  (T,  +  T2  )M(t )] ; 

h  x  Hinc  =  -h  x  [(T,  +  T2  )J(t )  +  (r,-%  +  rj;%  )M(t)] ; 

where  n  is  the  surface  normal  vector,  i/v=  1,2  the  medium  impedance,  and  L,  and  Tv  the  linear 
integral  operators  defined  in  reference  23.  It  should  be  noted  that  the  kernels  of  these  operators 
are  first  regularized  as  necessary  before  the  application  of  the  quadrature-based  scheme.  Then, 
the  interpolated  surface  profile  is  divided  into  equilength  segments  of  length  h  with  Nq 
quadrature  points  distributed  on  each  segment,  and  equation  1  is  transformed  into  a  matrix 
system  with  Gauss-Legendre  quadrature  as  the  underlying  quadrature  rule  of  the  LCN  method.  In 
the  present  formulation,  the  far- field  component  of  the  impedance  matrix  is  extracted  and  then  is 
characterized  by  a  standard  two-level  fast  multipole  method  (FMM)  of  0(N15)  complexity  in 
memory  and  central  processing  unit  (CPU)  load.  Intermediate-distant  interactions  are  accounted 
for  by  single-point  evaluations  of  the  kernels  in  equation  1.  Near  (or  local)  field  interactions  are 
provided  by  the  corrected  quadrature  weights  derived  from  testing  the  regularized  kernels  with 
Legendre  polynomials  of  orders  n  =  0,  1,  2,  . . ,  Nq  -  1.  The  overall  algorithm  achieves 
controllable  exponential  error  convergence  of  the  form  0(hNq ) ;  consequently,  fewer  number  of 
unknowns  are  needed  to  attain  the  solution — at  a  given  level  of  accuracy — as  compared  to 
conventional  low-order  method  of  moments  (MoM)  schemes. 

In  the  outlined  integral  equation  approach,  the  physical  surface  must  be  truncated  into  a  finite 
simulation  domain.  For  the  radar  scattering  problem,  mechanisms  must  be  build  into  the  solver  in 
order  to  sufficiently  suppress  the  surface  end  effects.  This  is  especially  critical  for  the  calculation 
of  the  signal  in  the  backscattering  directions,  where  the  surface  scattering  components  can  be 
much  weaker  than  the  edge  diffracted  waves.  In  this  study,  the  following  tapered  beam  (28)  is 
applied  as  the  excitation: 
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in  which  g  is  a  parameter  controlling  the  beamwidth  of  the  illumination.  As  equation  2  is 
essentially  only  a  closed-form  approximation  to  the  superposition  of  a  continuous  spectrum  of 
Gaussian-weighed  plane  waves,  it  does  not  satisfy  the  wave  equation  exactly — that  is,  it  is  not 
Maxwellian  in  nature.  In  fact,  for  fixed  g,  the  resultant  error  in  the  scattering  solution  grows  with 
decreasing  incident  angle.  It  has  been  shown  (28)  that  in  order  to  consistently  achieve  accurate 
solutions,  the  parameter  g  should  adhere  to  the  criterion 


g» 


1 


k,  sin  0. 


(5) 


A  more  substantive  definition  for  g  (or,  specifically,  its  minimum)  has  been  derived  by  Ye  and 
Jin  (29) — and  is  followed  by  the  present  work: 
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After  establishing  equation  2,  the  horizontal  surface  length  L  of  the  simulation  domain  should  be 
chosen  to  be  large  enough  to  ensure  the  surface  currents  near  the  edges  have  been  sufficiently 
attenuated.  It  is  seen  in  this  work  that  the  following  condition  should  be  met:  4 g  <  L  <  6 g,  where 
the  lower  limit  is  applicable  for  fast-varying  surface  profiles  and  the  upper  limit  for  gently 
varying  ones.  Because  of  the  requirements  imposed  on  the  tapered  beam,  it  is  apparent  that 
extremely  large  simulation  domains  are  needed  by  the  SIE  method  for  characterizing  grazing 
scattering  effects.  For  instance,  at  6inc  =  5°,  L  ~  14002.  Simulations  of  large  1-D  surfaces  at 
grazing  angles  have  been  carried  out  previously  (30,  31).  As  an  alternative  to  equation  2, 
excitation  in  the  form  of  an  integral  spectrum  of  tapered  plane  waves  has  been  proposed  (32,  id). 
Although  such  a  stimulus  is  fully  Maxwellian,  as  observed  in  this  work  and  others  (32,  33),  a 
large  simulation  domain  is  still  necessary  for  modeling  grazing  angle  phenomena,  and  therefore, 
its  implementation  does  not  provide  a  noticeable  computational  advantage  in  solving  equation  1 . 
A  partial  explanation  to  this  fact  can  be  reached  by  noting  that,  fundamentally,  the  surface 
supports  a  ground  wave  that  is  only  slightly  attenuated  by  the  roughness  even  though  the  rms 
slope  of  the  surface  variations  might  be  on  the  order  of  a  wavelength  or  more  (20).  Thus,  some 
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energy  propagates  out  to  the  ends  of  the  surface  (and  subsequently  gets  diffracted)  even  when  the 
ends  are  not  directly  illuminated  because  of  the  explicit  tapering.  A  complex  source  beam  also 
has  been  suggested  as  an  excitation  (22),  but  its  use  does  not  seem  to  provide  efficiency 
improvements  in  the  solver.  Instead  of  manipulating  the  incident  wave,  Oh  and  Sarabandi  (34) 
proposed  placing  resistive  sheets  at  the  ends  of  a  perfectly  conducting  surface  to  absorb  and 
minimize  the  induced  edge  currents;  however,  this  method  has  not  been  verified  for  dielectric 
surfaces  at  grazing  angles.  (The  resistive  loading  approach  has  been  improved  by  Zhao  and  West 
[35]  for  extension  to  surfaces  modeled  with  an  impedance  boundary  condition,  but  the  grazing 
angle  is  restricted  to  be  larger  than  20°.)  A  promising  technique  that  avoids  the  use  of  both 
excitation  and  resistive  tapering  has  been  implemented  by  Spiga  et  al.  The  so-called  grazing 
MoM  approach  (36)  applies  a  modified  integral  equation  formulation  to  a  locally  perturbed  half 
plane.  The  computational  domain  size  is  shown  to  be  independent  of  the  wave  incidence  angle. 


2.2  FDTD  Solver 

The  FDTD  method  directly  solves  the  differential  forms  of  Maxwell’s  equations  in  the  spatial 
and  temporal  domains.  Finite-difference  approximations  to  the  curl  expressions  for  the  electric 
and  magnetic  fields  generate  the  following  update  equations  (for  example,  for  vertical 
polarization,  and  assuming  isotropic,  non-dispersive,  non-permeable  media): 
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In  particular,  in  this  section,  the  2-D  solver  from  reference  37  is  used  as  the  basis  for  developing 
the  framework  needed  for  characterizing  scattering  from  the  dielectric  rough  surface.  In  the 
approach  here,  a  finite-length  rough  interface  between  two  media  is  embedded  in  an  infinite  half¬ 
space  environment  (i.e.,  two  dielectric  background  media  separated  by  a  flat  plane).  The  FDTD 
algorithm  implements  the  split  field  approach  (38),  in  which  the  computational  domain  is 
partitioned  into  total  and  scattered  field  regions  separated  by  a  connective  boundary  (figure  2).  In 
the  absence  of  any  scatterers  (i.e.,  only  the  half-space  background  is  present  in  the  scene),  the 
scattered  field  is  null  everywhere.  A  scatterer  is  defined  as  any  perturbation  to  the  half-space 
configuration  (such  as  a  roughness  of  the  interface)  and  generates  non-zero  scattered  field.  In  the 
case  of  a  rough  surface,  this  field  represents  the  incoherent  component  of  the  energy  scattered  by 
the  interface.  Note  that  the  incident  field  is  implemented  as  a  uniform  plane  wave  (fully 
Maxwellian),  with  the  sources  at  infinity.  The  incident  wavefront  is  not  tapered,  but  rather  a 
gradual  transition  from  the  ends  of  the  rough  surface  to  the  infinite  flat  ground  plane  is 
introduced  in  the  interface  profile.  As  such,  the  support  of  the  scattered  field  is  bounded  in  the 
numerical  domain,  and  the  method  is  akin  to  the  aforementioned  grazing  MoM  technique.  It  is 
important  to  emphasize  that  this  bounded  surface  perturbation  approach  is  different  from  those  of 
previous  FDTD  studies  on  rough  surface  scattering  (39-41)  in  that  only  the  perturbed 
(incoherent)  field  due  to  the  interface  roughness  is  calculated  here.  In  fact,  the  coherent 
component  of  the  field  reflected  by  the  surface  is  not  included  in  the  solution. 
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Figure  2.  Description  of  the  FDTD  computational  domain  with  a  bounded  surface  perturbation  for 
rough  surface  scattering  modeling. 
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One  major  issue  with  the  FDTD  method  is  the  so-called  numerical  dispersion  (38),  which  is  a 
result  of  the  finite-difference  approximation  of  the  field  differential  equations.  The  phase  errors 
introduced  by  this  effect  accumulate  as  the  wave  propagates  through  the  computational  domain 
and  can  become  significant  for  large  propagation  distances.  One  direct  consequence  is  that  in 
practice,  the  scattered  fields  produced  by  the  unperturbed  half-space  configuration  are  not 
exactly  null,  therefore  potentially  giving  rise  to  errors  that  may  be  comparable  in  magnitude  to 
the  incoherent  field  scattered  by  the  rough  surface.  In  order  to  minimize  the  numerical  dispersion 
artifacts,  the  unperturbed  solution  obtained  for  a  flat  interface  (or  the  numerical  dispersion 
error”)  is  coherently  subtracted  from  the  overall  solution  obtained  in  the  presence  of  the  rough 
surface. 

Throughout  this  study,  the  excitation  is  provided  by  an  unmodulated  fourth-order  Rayleigh  pulse 
(42).  Square  Yee  cells  are  used  to  discretize  the  entire  computational  domain;  thus,  the  interface 
profile  is  quantized  by  a  stair-stepped  approximation.  The  consequences  of  this  approximation 
on  rough  surface  modeling  are  discussed  in  section  2.3.  The  horizontal  extent  L  of  the  surface  is 
chosen  to  be  10  times  its  correlation  length — independent  of  the  incidence  angle.  To  emulate  an 
infinite  propagation  domain,  the  lattice  is  surrounded  by  a  perfectly  matched  layer  (PML)  on  all 
sides.  To  compute  the  far-field  scattering  pattern,  a  near-to-far  zone  transformation  is  performed 
with  the  scattered  field  sampled  at  points  situated  immediately  outside  the  connective  boundary 
(on  the  closed  Huygens  surface  in  figure  2). 

As  a  whole,  the  technique  outlined  previously  is  called  here  the  bounded  surface  perturbation 
with  coherent  error  subtraction  approach.  Its  region  of  validity  in  terms  of  the  surface  parameters 
and  incidence  angle  is  investigated  in  section  2.3.  As  opposed  to  SIE,  note  the  FDTD  approach 
requires  a  discretization  of  the  entire  space  enclosed  by  the  domain  boundary,  and  additionally 
demands  the  sampling  to  be  compatible  with  the  highest  frequency  and  a  domain  size  compatible 
with  the  lowest  frequency.  Nevertheless,  the  low  implementation  complexity  and  the  ability  to 
easily  model  ground  (and  target)  inhomogeneities — coupled  with  the  fact  that  the  required  size 
of  the  simulation  domain  is  independent  of  the  incident  angle — make  the  FDTD  algorithm  a 
more  attractive  option  for  characterizing  the  grazing  scattering  behaviors  of  terrain  scenes  as 
needed  in  this  work. 

2.3  Comparison  of  SIE  and  FDTD  Results 

For  communication  systems  propagation  scenarios,  the  accuracy  of  the  SIE  solver  has  been 
verified  in  reference  23.  For  grazing-angle  scattering  problems,  convergence  studies  (in  terms  of 
g  and  L)  have  been  carried  out  within  this  work  to  validate  the  solver’s  robustness.  Subsequently, 
as  described  in  this  section,  efforts  are  taken  to  benchmark  the  accuracy  of  the  proposed  FDTD 
algorithm  using  the  SIE  solutions. 

In  both  the  SIE  and  FDTD  solvers,  the  ground  is  modeled  as  a  homogeneous  medium  with  an 
effective  relative  dielectric  constant  sr  and  conductivity  cr(/.  In  general,  the  surface  statistics  of  the 
ground  are  described  by  two  functions:  the  probability  density  function  of  the  height  variations 
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and  the  surface  autocorrelation  function.  Within  the  scope  of  this  study,  a  zero-mean  surface 
profile  s(x)  obeying  Gaussian  statistics  {hrms,  lc} — where  hrms  is  the  rms  height  and  lc  the 
correlation  length — is  assumed  and  generated  from  its  randomized  power  spectral  density 
function  by  following  the  procedure  prescribed  in  reference  28.  For  the  S1E  solver,  as  it  is 
impractical  to  use  s(x)  directly  because  of  the  necessity  for  rapid  generation  and  inter-translation 
of  surface  height,  length,  and  slope  parameters,  an  approximate  form  for  s(x),  which  facilitates 
efficient  computation  of  these  routines,  must  be  employed.  Therefore,  an  approximate  s(x )  is 
reconstructed  from  a  sampled  set  of  s(x)  using  the  approximate  prolate  spheroidal  wave  function 
(APSWF)  (43)  as  the  basis.  Although  only  surfaces  with  Gaussian  statistics  are  considered  in  this 
work,  the  tools  described  are  flexible  enough  to  allow  the  generation  and  analysis  of  arbitrary 
band-limited  roughness  profiles. 

In  interpreting  the  simulated  responses,  the  radar  parameter  of  interest  is  the  incoherent 
scattering  coefficient  of  the  surface,  which  is  approximated  here  as 

(7°  =  lim— - , - (10) 

'-*0  I u.  I  L.„ 

for  the  FDTD  solver,  Lmu,  the  effective  length  of  the  surface  illuminated  by  the  incident  field,  is 
the  entire  length  of  the  surface  perturbation,  whereas,  for  the  SIE  solver,  Lmu  =  g(n!2)  .  The 
ensemble  average  in  reference  10  is  calculated  using  60  surface  realizations  for  each  set  of 
surface  parameters  {hrms,  lc) .  Monte  Carlo  simulations  are  carried  out  at  500  MHz  for  hrms  =  0. 12 
to  12,  lc  =  0. 12  to  12  in  steps  of  0. 12.  As  most  realistic  terrain  surfaces  have  an  rms  height  that  is 
less  than  the  correlation  length,  to  reduce  the  number  of  simulation  runs,  only  surface  sets  with 
hrms  <  lc  are  considered.  Comparisons  of  the  FDTD-  and  SIE-calculated  backscattering 
coefficient  for  6inc  =  5°,  10°,  15°,  and  45°  are  shown  in  figures  3  and  4,  for  vertical  and 
horizontal  polarization,  respectively.  The  error  in  the  FDTD  results  is  defined  as 
A  =  <7°dB  FDTD  -  (J°B  SIE  .  Note  that  over  the  range  of  surface  parameters  of  interest,  very  good 

agreement  is  observed  between  the  two  solvers  for  both  polarizations,  except  for  the  region  of 
small  rms  height  and  large  correlation  length  (on  the  lower  right  comer  of  the  figures — the 
-inaccurate”  region)  where  the  FDTD  solution — at  first  glance — seems  to  deviate  from  the  SIE 
solution.  Outside  this  region,  for  vertical  polarization,  the  average  errors  ( |  A| )  are  1.4  dB,  1 .2  dB, 

1.2  dB,  and  0.8  dB  for  6inc  =  5°,  10°,  15°,  and  45°,  respectively.  Similarly,  for  horizontal 
polarization,  the  average  errors  are  1.4  dB,  1.3  dB,  1.1  dB,  and  1.0  dB.  Although  the 
backscattering  coefficient  is  the  primary  quantity  of  relevance  in  this  study,  for  a  more  detailed 
comparison  of  the  two  solvers,  the  bistatic  scattering  patterns  of  the  rough  surface  are  analyzed 
in  figures  5  and  6.  Surface  parameter  sets  in  a  region  where  FDTD  is  well  behaved  and  in  one 
where  it  deviates  from  the  SIE  solution  are  selected  for  comparison.  It  is  seen  that  even  for  cases 
where  erroneous  results  are  observed  in  the  backscattering  directions,  the  FDTD  routine  still 
retains  an  accurate  solution  in  the  forward  scattering  directions — this  stems  from  the  fact  that  the 
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relative  error  is  larger  in  the  directions  of  weak  scattering  response,  but  the  absolute  error  is 
likely  the  same  throughout  the  range  of  observation  angles. 


Figure  3.  Errors  (A  in  dB)  in  the  FDTD-calculated  vertical-vertical  (vv)-polarized  backscattering 

coefficient  (  G°^  ) — as  validated  with  the  SIE  solution — at  various  incidence  angles  for  ground 

profiles  characterized  by  Gaussian  statistics  {hrms,  lc };  sr  =  5.56,  od  =  10  mS/m.  Annotations  at 
data  points  are  shown,  e.g.,  at  6inc  =  5°,  hrms  =  0.32,  lc  =  0.72,  A  =  0.40  dB. 
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Figure  4.  Errors  (A  in  dB)  in  the  FDTD-calculated  horizontal-horizontal  (M)-polarized  backscattering 
coefficient  ( <J°hh  ) — as  validated  with  the  SIE  solution — at  various  incidence  angles  for  ground 
profiles  characterized  by  Gaussian  statistics  {hrms,  lc};  sr  =  5.56,  ad  =  10  mS/m. 
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Figure  5.  Comparison  of  vv-polarized  bistatic  scattering  patterns  as  derived  by  SIE  (solid  lines)  and  FDTD  (dash 
lines)  at  various  incident  angles  for  two  sets  of  surface  parameters:  {0. I  A,  \A)  and  {0.3/1,  U};  sr  =  5.56, 
od  =  1 0  mS/m. 
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Figure  6.  Comparison  of  M-polarized  bistatic  scattering  patterns  as  derived  by  SIE  (solid  lines)  and  FDTD  (dash 
lines)  at  various  incident  angles  for  two  sets  of  surface  parameters:  {0. \1,  \1}  and  {0.31,  11};  er  =  5.56, 
od=  10  mS/m. 


For  all  the  simulations  carried  out  in  figures  3-6,  the  cell  size  is  uniformly  chosen  to  be  0.01k — 
though  it  is  seen  that  a  much  coarser  cell  size  can  be  used  for  surfaces  with  parameters  located 
away  from  the  inaccurate”  region.  Figure  7  demonstrates  that  the  aforementioned  inaccuracies 
in  the  backscattering  direction  of  the  FDTD  solution — for  surfaces  with  small  rms  slope  and 
large  correlation  length — can  be  reduced  by  further  refining  the  FDTD  grid.  A  smaller  cell  size, 
of  course,  means  increased  computational  burden.  It  should  be  mentioned  that  increasing  the 
length  of  the  simulation  domain  (e.g.,  from  10/c  to  20 lc)  does  not  significantly  improve  the 
solution  accuracy.  Interestingly,  the  extremely  weak  backscattering  responses  generated  by 
surfaces  with  small  rms  slopes  also  pose  a  challenge  for  the  SIE  solver,  as  the  length  of  the 
simulation  domain  must  be  made  large  enough  to  guarantee  the  suppression  of  the  edge 
diffraction  effects.  For  instance,  as  seen  in  this  study,  for  a  surface  with  parameters  {0.3k,  lk},  a 
surface  length  of  4 g  is  adequate;  whereas,  for  a  surface  with  parameters  {0.1k,  lk},  a  length  of 
6g  is  needed. 
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Figure  7.  Convergence  of  FDTD  solution  as  function  of  cell  size  for  surface  with 
parameters  {0. 1A,  IA};  sr  =  5.56,  ad=  10  mS/m;  0inc  =  5°. 
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3.  Extension  of  FDTD  Simulation  of  Rough  Surface  Scattering  to  3-D 


Though  the  solution  to  the  2-D  scattering  problem  (for  the  1-D  surface)  can  convey  important 
physical  insights,  only  the  full-wave  vectorized  solution  to  the  3-D  scattering  problem  (for  the 
2-D  surface)  can  provide  a  complete  understanding  of  the  scattering  and  propagation 
mechanisms,  including  cross-polarization  interactions.  Integral  equation-based  3-D  fast  solvers 
devoted  to  analyzing  rough  surface  effects  have  been  proposed  {44-47).  For  instance,  the 
SDFMM  {47)  achieves  0{N)  complexity  for  both  CPU  time  and  memory  requirements,  as 
opposed  to  an  unaccelerated  scheme,  which  has  0(N  )  complexity.  As  reported  in  literature, 
these  fast  algorithms  have  been  applied  for  simulating  only  relatively  small  surfaces  with  non¬ 
grazing  angle  excitations.  However,  it  should  be  mentioned  that  surfaces  with  dimensions  up  to 
2561  x  321  have  been  simulated  by  Torrungrueng  and  Johnson  {48)  in  analyzing  low-grazing 
angle  effects.  High-order  routines  and  phase-extraction  (23)  can  be  used  to  explicitly  reduce  the 
number  of  unknowns  in  the  integral  equation  method;  however,  no  works  have  been  done  yet 
that  take  advantage  of  these  techniques  for  the  3-D  grazing  scattering  problem. 

As  complements  to  SIE  routines,  FDTD  algorithms  for  studying  scattering  from  2-D  rough 
surfaces  have  also  been  developed.  In  reference  49,  the  total  field/scattered  field  (TF/SF) 
approach  is  used  for  zoning  a  computational  lattice,  which  is,  in  turn,  truncated  by  averaging 
boundary  conditions  instead  of  a  PML.  In  the  scattered  field  calculation  stage,  a  windowing 
function  is  applied  to  the  near-to-far-field  transformation  surface  to  reduce  the  boundary 
reflection  effects.  In  reference  50,  exploiting  periodic  boundary  conditions,  the  2-D  surface  is 
modeled  as  a  periodic  structure  with  one  period  of  the  surface  extracted  for  FDTD  computations. 
As  the  update  equations  in  FDTD  are  inherently  parallel  in  nature,  Message-Passing  Interface 
(MPI)-based  routines  also  have  been  implemented  to  facilitate  the  simulation  of  large  domains. 
The  TF/SF  algorithm  in  reference  5 1 — developed  specifically  for  studying  rough  surface 
scattering — adopts  the  uniaxial  PML  for  terminating  the  computational  domain  and — similar  to 
what  is  done  in  SIE  routines — employs  a  Gaussian  tapered  incident  wave  to  suppress  truncation 
effects.  Unfortunately,  the  focus  of  the  aforementioned  routines  has  been  on  non-grazing  angle 
illuminations  and  the  effectiveness  of  the  various  truncation  boundary  conditions  and  incident 
wave  treatments  for  modeling  the  low-grazing  angle  scattering  problem  has  not  been  studied. 

In  this  work,  the  modeling  of  rough  surface  scattering  is  extended  to  the  3-D  configuration  by 
employing  the  AFDTD  code.  AFDTD  (52)  is  a  software  package  developed  at  ARL  for  radar 
signature  calculation  and  is  based  on  a  3-D  implementation  of  the  FDTD  algorithm.  The  code  is 
fully  parallelized  and  runs  on  large  distributed  computer  systems  using  the  MPI  framework.  The 
computational  domain  is  decomposed  into  rectangular  subdomains  and  the  FDTD  equations  are 
solved  separately  for  each  subdomain  within  one  MPI  process.  Only  the  electric  and  magnetic 
field  samples  in  the  boundary  layer  between  two  adjacent  subdomains  need  to  be  exchanged 
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between  the  respective  processes  (55).  As  a  consequence,  the  AFDTD  code  is  highly  scalable, 
even  when  the  number  of  MPI  processes  is  in  the  hundreds.  The  rough  surface  scattering 
implementation  for  the  3-D  AFDTD  code  follows  closely  the  bounded  surface  perturbation 
approach  presented  in  section  2  for  the  2-D  version.  Random  rough  surfaces  with  isotropic 
statistical  properties  are  examined  in  this  section.  Therefore,  on  average,  the  scattered  power  in 
the  plane  of  incidence  is  independent  of  the  azimuth  incidence  angle. 

The  expected  range  of  validity  (in  terms  of  surface  parameters  and  cell  size)  of  the  FDTD 
method  has  been  examined  in  section  2.  Further  validation  is  performed  by  comparing  3-D 
Monte  Carlo  scattering  results  to  measured  data  reported  in  reference  9.  Two  types  of  ground 
surface  are  considered:  for  figure  8,  hrms  =  4  mm,  lc  =  8.4  cm,  sr  =  15.57,  and  Od  =  0.31  S/m;  and 
for  figure  9,  hrms  =  3  cm,  lc  =  8.8  cm,  er  =  8.92,  and  oj  =  0.19  S/m.  The  figures  show  there  is  good 
agreement  between  the  simulated  and  measured  backscattered  responses  (at  1.5  GHz)  for 
20°  <  6inc  <  80°  for  both  co-polarization  components.  Note  the  measured  response  includes  a 
coherent  component  (which  accounts  for  the  strong  peak  observed  in  figure  8  near  normal 
incidence)  whereas  the  FDTD  response  does  not.  Also,  while  the  FDTD-simulated  surfaces  are 
purely  Gaussian,  the  measured  profiles  are  of  a  more  complex  spectrum  and  can  be  fitted  only 
approximately  by  Gaussian  statistics — this  perhaps  explains  the  discrepancy  seen  near  6inc  =  20° 
for  G°hh  in  figure  8.  Unfortunately,  no  empirical  data  can  be  found  in  literature  at  grazing  angles 

less  than  20°  over  the  interested  frequency  range.  The  ensemble  averages  of  the  results  from  30 
surface  realizations  are  shown  in  these  figures:  each  surface  realization  has  dimensions 
10/c  x  10/C;  solution  convergence  is  observed  with  a  FDTD  grid  size  of  1  mm;  the  simulation  of 
one  realization  at  each  incidence  angle  takes  approximately  3  h  on  four  2.8-GHz  Intel  Nehalem 
processors,  with  a  total  memory  load  of  6.6  GB. 
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Figure  8.  Comparison  of  FDTD-simulated  3-D  backscattering  results  with  measured  data  from 
reference  9:  hrms  =  4  mm,  lc  =  8.4  cm,  sr  =  15.57,  and  od  =  0.31  S/m.  At  angles  close  to 
normal  (0inc  —>  90°),  the  strong  measured  response  is  due  to  the  inclusion  of  a  coherent 
surface  reflection  component;  this  component  is  not  modeled  in  the  simulations. 
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Figure  9.  Comparison  of  FDTD-simulated  3-D  backscattering  results  with  measured  data 
from  reference  9:  hrms  =  3  cm,  lc  =  8.8  cm,  er  =  8.92,  and  ad  =  0. 19  S/m. 


The  developed  numerical  tool  is  invaluable  for  studying  the  angular  and  spectral  behaviors  of 
ground  clutter.  Displayed  in  figure  10  is  the  backscattering  coefficient — over  the  frequency  band 
500-2500  MHz,  for  5  °  <  f)mc  <15  0 — of  a  typical  surface  with  hnns  =  1  cm,  lc  =  7  cm,  er  =  5.56, 
and  Od  =  10  mS/m.  It  is  observed  that,  in  general,  the  backscattering  power  increases  with  the 
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grazing  angle  and  cf'v  >  <y°hh  >  <j°hv  ( <y°vh  =  <y°hv  by  reciprocity).  Maximum  surface  scattering 

appears  to  be  concentrated  within  900-1500  MHz  in  this  example.  (The  non-smooth  nature  of 
the  intensity  variations  within  the  figures  is  a  numerical  artifact  of  the  Monte  Carlo  simulations.) 
Qualitatively,  the  detectability  of  a  specific  target  can  be  determined  by  simply  comparing  its 
scattering  coefficient  to  that  of  the  ground  surface.  A  more  complete  consideration  requires  the 
intensity  distribution  of  the  ground  response  (in  the  image  domain)  to  be  studied.  For  the  band  of 
frequencies  examined  here,  figure  10  shows  that  the  amplitude  of  the  surface  response  is 
increasing  with  frequency  at  the  lower  end  of  the  band  but  decreasing  with  frequency  at  the 
higher  end.  This  behavior  is  characteristic  to  Gaussian  surfaces  and  consistent  with  results 
presented  in  a  previous  experimental  study  (54).  Note  that  some  real  surface  profiles  may  be  of  a 
more  complex  spectrum  and  can  be  only  approximated  by  Gaussian  statistics.  In  fact,  the 
autocorrelation  function  of  real  natural  terrain  surfaces  can  be  approximated  to  be  Gaussian, 
exponential,  or  a  combination  of  the  two  (9).  The  Gaussian  surface  is  expected  to  be  an  adequate 
model  at  low  frequencies  where  the  fast  surface  variations  inherent  in  an  exponential  spectrum 
are  not  sensed  by  the  wave. 
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Figure  10.  Angular  and  spectral  distribution  of  backscattering  from  a  rough  surface 

with  surface  parameters  hrms  =  1  cm,  lc  =  7  cm,  er  =  5.56,  and  ad=  10  mS/m. 
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4.  Emulation  of  Forward-Looking  Radar  Imaging  in  Presence  of  Ground 
Clutter 


Having  studied  the  accuracy  and  region  of  validity  of  the  proposed  FDTD  solver,  large-scale 
full-wave  modeling  of  multistatic  target  imaging  in  a  rough  ground  environment  is  described  in 
this  section.  An  investigation  of  particular  interest  is  the  analysis  of  the  statistics  of  the  ground 
response  in  the  images. 

4.1  Emulation  and  Imaging  Framework 

First,  the  parallelized  3-D  FDTD  algorithm  is  applied  to  simulate  composite  scattering  from 
targets  in  a  rough  ground  environment.  Then,  a  coherent  field  integration  technique  is  employed 
to  process  the  scattered  signals  to  obtain  an  image  of  the  illuminated  scene.  According  to  the 
time-reversal,  or  phase-conjugation,  method  (55,  56),  since  the  received  scattered  field  can  be 
written  as 

Es(fR,fT,a)  =  T((o)a{o))G{fs,fT,o})G(fR,fs,(o),  (11) 

an  approximate  objective  (image)  function  can  be  formed  with  the  following  expression: 

co  R  T 

where  the  asterisk  notation  denotes  the  phase  conjugation  operation;  G(r,r,a> )  is  the  Green’s 

function  of  the  environment  ( r  is  the  position  of  the  observation  point  and  r  is  the  position  of 
the  source  point);  rT ,  rR,  and  r  are  the  locations  of  the  transmitter,  receiver,  and  scatterer, 

respectively;  and  T((o)  and  a(  co)  are  the  spectra  of  the  transmitted  waveform  and  target 

response.  The  angular  and  range  resolution  of  the  image  are  determined  by  the  sensor  aperture 
and  system  bandwidth,  respectively,  of  the  coherent  summation  performed  in  equation  12.  For 
sensing  in  the  presence  of  a  randomly  varying  ground  interface,  as  the  exact  propagation  Green’s 
function  is  not  known,  G(r,r',co)  is  approximated  by  the  half-space  Green’s  function.  As  such, 

the  target  image  is  corrupted  by  ground  clutter. 

The  effects  of  a  rough  ground  surface  on  focusing  are  demonstrated  with  the  simulation  results 
displayed  in  figures  11-13,  which  show  images  of  a  9  m  x  19  m  area  populated  with  targets  in 
the  form  of  landmines  (either  anti-tank  or  anti-personnel,  metallic  or  plastic,  buried  or  on- 
surface)  and  155-mm  shells  (metallic,  buried).  Specifically,  the  targets — the  shapes,  dimensions, 
orientations,  and  locations  of  which  are  indicated  by  the  white  outlines  in  image  (a) — include  the 
following: 

1 .  Buried  metallic  anti-personnel  landmine; 
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2.  Buried  plastic  (er  =  3. 1,  Od  =  2  mS/m)  anti-personnel  landmine; 

3.  On-surface  plastic  (er  =  3.1,  rrt/  =  2  mS/m)  anti-personnel  landmine; 

4.  Buried  metallic  155-mm  shell; 

5.  Buried  metallic  anti-tank  landmine; 

6.  On-surface  metallic  anti-tank  landmine; 

7.  Buried  metallic  155-mm  shell; 

8.  Buried  metallic  155-mm  shell; 

9.  On-surface  metallic  anti-personnel  landmine; 

10.  Buried  plastic  (er  =  3. 1,  Od  =  2  mS/m)  anti-tank  landmine; 

11.  On-surface  plastic  {er  =  3.1,  rrt/  =  2  mS/m)  anti-tank  landmine. 
Buried  targets  are  positioned  at  3  cm  below  the  surface. 
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Figure  11.  Simulated  vv-polarized  images  for  9  m  x  19  m  scene:  (a)  Ground  with  flat  surface;  (b)  ground  with 
randomly  rough  surface,  hrm  =  1.2  cm,  lc  =  14.93  cm;  and  (c)  ground  with  randomly  rough  surface, 
hrms  =1.6  cm,  lc  =  14.93  cm.  Ground  electrical  properties:  sr=  8,  ad  =  10  mS/m. 
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Figure  12.  Simulated  M-polarized  images  for  9  m  x  19  m  scene:  (a)  Ground  with  flat  surface; 

(b)  ground  with  randomly  rough  surface,  hrm  =  1.2  cm,  lc  =  14.93  cm;  and  (c)  ground  with 
randomly  rough  surface,  hrms  =1.6  cm,  lc  =  14.93  cm.  Ground  electrical  properties:  sr  =  8, 
od  =  10  mS/m. 
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Figure  13.  Simulated  Av-polarized  images  for  9  m  x  19  m  scene:  (a)  Ground  with  flat  surface; 

(b)  ground  with  randomly  rough  surface,  hrm  =  1.2  cm,  lc  =  14.93  cm;  and  (c)  ground  with 
randomly  rough  surface,  hrms  =1.6  cm,  lc  =  14.93  cm.  Ground  electrical  properties:  er  =  8, 
=10  mS/m. 


At  this  point,  some  important  differences  between  the  configuration  of  the  SIRE  radar  and  the 
setup  of  the  simulation  scenario  should  be  mentioned.  The  SIRE  radar  transceiver  consists  of  an 
array  of  equally  spaced  receivers  placed  across  a  linear  aperture  and  two  transmitters  situated  at 
the  ends  of  the  aperture.  The  entire  assembly  is  elevated  at  2-m  height  (atop  a  vehicle)  and 
covers  downrange  between  8  and  21m  while  the  platform  is  on  the  move.  (To  be  specific,  the 
radar  height  can  be  varied  from  2  to  3  m;  in  field  experiments,  the  downrange  coverage  is  from 
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8  to  33  m,  although  only  the  captured  scattered  fields  corresponding  to  a  distance  of  8  to  21  m 
are  used  for  image  formation.)  Since  AFDTD  only  implements  incidence  and  scattering  in  the 
far-field,  plane-wave  incidence  from  two  azimuth  directions  and  five  elevation  angles  that  span 
the  same  solid  angle  as  the  SIRE  transmitters  is  considered.  For  each  elevation  angle,  the  bistatic 
scattered  fields  are  computed  at  equally  spaced  azimuth  angles  (in  1°  increments)  that  provide 
similar  azimuth  coverage  as  the  SIRE  receivers.  At  the  limits  of  the  SIRE  radar  range  coverage, 
the  elevation  angles  and  angular  apertures  correspond  to  6inc  =  5°,  A(f>  =  32°  and  0inc  =  15°, 

A</>  =  11°,  respectively.  The  FDTD  simulations  were  performed  at  the  ARL  Distributed 
Supercomputing  Resource  Center  (DSRC)  in  Aberdeen,  MD,  on  a  SGI  Altix  ICE  8200  system. 
The  discretization  of  the  entire  9  m  x  19  m  scene  involved  1.9  billion  cells,  and  each  simulation 
run  at  one  incidence  angle  required  256  processors  and  approximately  2000  CPU  hours. 

4.2  Results  and  Discussions 

The  co-  and  cross-polarization  images  generated  from  the  full-wave  simulation  data  are 
displayed  in  figures  11-13,  where  three  ground  roughness  scales  are  examined:  (a)  flat  ground; 
(b)  hrms  =1.2  cm,  lc  =  14.93  cm;  and  (c)  hrms  =1.6  cm,  lc  =  14.93  cm.  The  response  from  the 
larger  targets  (anti-tank  landmines)  in  the  backscattering  directions  is  primarily  due  to  edge 
diffractions  from  the  front  and  back  rims  of  the  top  of  the  cylindrical  structure.  The  scattering 
centers  corresponding  to  these  two  edges  are  resolved  in  the  images — for  both  on-surface  and 
buried  cases.  For  the  smaller  targets  (anti-personnel  landmines),  the  two  scattering  edges  tend  to 
merge  into  one,  since  the  imaging  configuration  as  chosen  here  is  unable  to  provide  sufficient 
resolution  to  separate  the  two  phase  centers.  In  the  presence  of  a  rough  ground  interface,  these 
diffraction  effects  are  less  apparent  in  the  images  for  the  buried  targets  due  to  interference  from 
surface  scattering.  Buried  targets  are  generally  more  difficult  to  discern  from  the  background  as 
compared  to  on-surface  targets.  For  instance,  the  co-polarized  intensity  of  the  buried  metallic 
anti-tank  landmine  is  about  15  to  20  dB  below  that  of  the  on-surface  configuration — for  a  ground 
with  either  a  flat  or  rough  surface.  The  difference  for  the  cross-polarized  return  is  less,  at  8  to 
10  dB.  A  buried  plastic  target  (either  large  or  small  landmine)  is  especially  hard  to  detect  owing 
to  the  limited  dielectric  contrast  between  the  target  and  the  soil  background.  The  orientation  of  a 
target  can  also  have  a  significant  impact  on  its  radar  return:  maximum  response  from  the 
155-mm  shell  is  observed  when  it  is  oriented  parallel  to  the  imaging  aperture.  A  much  weaker 
response  is  noted  when  it  is  oriented  perpendicular  or  obliquely  (at  45°)  to  the  aperture.  These 
observations  are  consistent  with  the  analysis  previously  reported  by  reference  57. 

Comparing  all  target  responses  across  the  three  polarizations,  on  average,  vv  provides  the 
strongest  response,  and  it  exceeds  hh  by  about  10  dB  (except  for  the  large  on-surface  landmines). 
The  hh  polarization,  in  turn,  exceeds  the  cross-polarization  response  by  about  7  dB.  For  the  large 
on-surface  landmines  (either  metallic  or  plastic),  the  two  co-polarization  components  are 
approximately  equal  in  intensity  due  to  grazing  angle  incidence.  (At  steeper  incidence  angles,  the 
hh  response  is  generally  larger  as  a  result  of  enhanced  ground  bounce  contribution).  However, 
when  these  targets  are  buried,  vv  again  seems  to  give  the  stronger  return.  As  a  vertically 
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polarized  signal  experiences  less  transmission  loss  at  the  air-ground  interface  and  suffers  less 
propagation  loss  in  general,  the  above  observation  supports  the  notion  that  vv  is  the  more 
desirable  configuration  for  detecting  targets  either  on  top  of  or  buried  under  a  flat  ground.  In  the 
presence  of  ground  surface  roughness,  detection  performance  must  be  evaluated  with  a 
consideration  of  the  target-to-background  clutter  ratio.  As  such,  since  the  surface  response  is 
strongest  at  vv,  vertical  polarization  does  not  seem  to  provide  an  advantage — as  seen  in  figures 
11c,  12c,  and  13c,  where  certain  buried  targets  (i.e.,  targets  4  and  5)  are  found  to  stand  out 
against  the  background  more  in  hh  and  horizontal-vertical  (hv)  than  in  vv. 

It  is  interesting  to  also  note  the  -butterfly-shaped”  formation  generated  by  each  point  target  in 
the  /n’-polarized  images:  two  pairs  are  observed  for  the  anti-tank  landmines  and  one  pair  for 
other  smaller  targets.  An  explanation  for  the  appearance  of  this  feature  can  be  derived  by 
realizing  that — for  a  vertically  polarized  (or  ^-polarized)  incident  wave — although  the  induced 
currents  on  the  targets  are  flowing  mostly  in  both  the  vertical  (parallel  to  z-axis)  and  downrange 
(parallel  to  x-axis)  directions,  only  the  downrange-oriented  currents  can  generate  any  cross- 
polarized  (Ef)  field  component.  As  a  result,  since  the  component  of  the  scattered  field  from 

these  downrange  currents  is  anti-symmetric  (that  is,  equal  in  amplitude,  opposite  in  sign)  about 
the  direction  of  the  incidence  wave,  integration  of  fields  due  to  excitation  from  two  incidence 
angles  (i j>\  and  <j> 2)  over  the  same  angular  span  A<f>  leads  to  the  formation  of  a  null  at  the  center  of 
the  target  image  pattern.  For  point-like,  discrete  targets,  although  symmetry  in  the  target 
structure  helps  in  establishing  a  stronger  response  for  this  — btierfly-shaped”  pattern,  it  is  not  a 
necessary  condition:  the  same  pattern  is  observed  for  the  155-mm  shell  (targets  4,  7,  and  8)  as 
well.  This  distinctive  feature  mentioned  could  be  exploited  for  target  identification  in  detection 
algorithms,  but  one  should  note  that  a  very  rough  ground  surface  could  exhibit  many  peaks 
acting  as  strong  point  scatterers  and,  in  turn,  generating  the  same  type  of  feature. 

The  images  of  figures  1 1-13  essentially  map  the  spatial  variation  of  the  average  radar  cross 
section  (RCS)  of  the  scene,  where  the  averaging  operation,  equation  12,  is  being  carried  over  a 
range  of  frequencies  and  observation  angles  in  elevation  and  azimuth.  An  overlay  of  the  images 
on  the  physical  surface  profile  suggests  that  the  brightest  surface  returns  are  from  the  surface 
peaks — or,  more  specifically,  from  the  faces  of  the  surface  peaks  that  are  oriented  toward  the 
radar.  This  is  different  from  radar  observations  at  larger  incidence  angles,  for  which  both  the 
peaks  and  troughs  of  the  surface  could  produce  strong  returns.  As  the  roughness  scale  increases, 
at  grazing  angles,  it  is  also  expected  that  the  shadowing  effect  would  play  a  significant  role  in 
determining  the  total  surface  return.  Curiously,  in  reviewing  the  rough  ground  image  data,  the  vv 
response  is  seen  to  be  much  greater  than  both  the  hh  and  the  cross-polarized  responses  as 
expected,  but  the  hh  response  itself  is  only  slightly  greater  than  the  cross-polarized  response.  For 
example,  according  to  figures  1  lb,  12b,  and  13b,  the  mean  ground  clutter,  cf(dB ,  is  calculated  to 

be  -55. 12  dB  for  vv,  -71.50  dB  for  hh,  and  -74.60  dB  for  hv,  while  a  study  of  the  backscattering 
coefficients  (within  the  same  frequency  band)  for  the  same  ground  properties  reveals  that  the  hh 
response  is  10  to  19  dB  higher  than  the  hv  response.  (The  large  difference  between  the 
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amplitudes  of  the  vv  and  hh  responses  is  a  result  of  observation  at  near-grazing  angles  and  the 
use  of  Gaussian  statistics  for  the  ground  surface.  In  view  of  the  lack  of  experimental  data  at  these 
grazing  angles,  it  is  unclear  whether  this  large  difference  is  also  evident  in  actual  measurement 
results  for  real  surfaces  of  a  more  complex  spectrum.)  The  discrepancy  between  the  imaging  and 
backscattering  results  can  be  resolved  by  analyzing  the  bistatic  scattering  pattern  of  the  rough 
surface  and  noting  that  the  cross-polarized  response  can  become  greater  than  the  hh  response 
away  from  the  incidence  direction  (f)inc.  In  fact,  at  low  elevation  angles  and  over  the  range  of 
azimuth  angles  (A ,<j>  <  90°)  of  interest  in  this  study,  the  hv  scattering  pattern  has  a  minimum  at 
<j)inC  but  quickly  increases  in  intensity  away  from  that  direction.  On  the  other  hand,  the  hh 
response  tends  to  remain  constant  or  even  decrease  over  the  same  range  of  observation  angles. 
Consequently,  the  average  response  taken  over  A <f>  may  lead  to  a  higher  cross-polarized  image 
response  relative  to  hh  than  as  indicated  by  backscattering  analysis  alone.  For  grazing-angle 
observations,  it  is  noted  that  these  behaviors  of  the  hh  and  cross-polarized  responses  (in  azimuth) 
are  equally  applicable  to  an  on-surface  or  shallow-buried  point-like,  discrete  target. 

The  vertical-horizontal  (v/z)-polarized  images  for  the  scene  examined  in  figures  1 1-13  are  not 
included  here  but  are  seen  to  be  very  similar  to  the  hv  case.  It  should  be  mentioned  that  given  the 
bistatic  sensing  geometry  used  in  this  work,  the  reciprocity  principle  cannot  be  exploited  for 
direct  inference  of  the  vh  response  from  the  hv  one,  or  vice  versa. 

To  quantify  target  detection  performance,  the  distribution  of  the  ground  response  needs  to  be 
estimated.  Once  the  image  clutter  distribution  is  known  and  an  acceptable  false-alarm  rate 
chosen,  a  threshold  value  for  the  detectability  of  a  target  of  a  given  RCS  can  be  determined.  For 
a  statistically  homogeneous  terrain  scene,  in  the  absence  of  multipath  and  shadowing  effects,  the 
scattered  field  response  (i.e.,  equivalently,  the  amplitude  of  the  scattering  matrix  elements,  |AP?|, 
or  radar  voltage  response)  is  commonly  found  to  conform  to  Rayleigh  statistics  (58).  In  other 
words,  the  real  and  imaginary  parts  of  Spq  are  observed  to  be  uncorrelated  zero-mean  Gaussian 
random  variables  with  equal  variances.  Other  distributions  also  have  been  proposed  to  account 
for  terrain  scenes  exhibiting  more  complex  and  non-Gaussian  scattering  phenomena.  Examples 
include  the  lognormal,  Weibull,  and  ^-distributions.  A  study  of  the  surface  clutter  in  the 
simulated  images  derived  in  this  work  indicates  that  while  a  Rayleigh  fading  model  is  adequate 
for  surfaces  with  small  rms  slopes,  a  ^-distribution  is  seen  to  be  more  appropriate  at  larger  rms 
slopes.  As  the  ^-distribution — which  is  a  mixture  of  a  Rayleigh  and  Gamma  distribution — can 
be  reduced  to  Rayleigh  as  a  special  case  (59),  the  ^-distribution  is  used  henceforth  as  the  basis  to 
model  the  rough  ground  response. 

In  the  image  domain,  the  probability  density  of  the  pixel  intensities  (i.e.,  the  scattered  power 
response  <jcdB ,  expressed  in  decibel)  is  found  to  be  a  very  good  match  to  the  function 
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where  -1  <  v  <  oo;  £  >  0;  s  =  j An ;  T  (•)  is  the  Gamma  function;  and  Kv  (•)  is  the  modified 

Bessel  function  of  the  second  kind  of  order  v.  The  parameter  v  controls  the  shape  of  the  density 
function,  with  £  as  a  scaling  factor.  It  is  important  to  note  that  <jcdB  itself  is  not  ^-distributed, 

though  expression  13  is  derived  here  by  assuming  a  ^-distribution  for  the  scattered  field 
amplitude.  Also,  note  for  large  argument  v,  it  can  be  shown  that 

p(<?dB  I  v large,  £)  *——•  —  £?  CT  ,  (14) 
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Expression  14  is  the  distribution  corresponding  to  Rayleigh-based  fading  statistics  with  mean 
power  (jc  (in  linear  scale).  As  the  /^-distribution  model  was  originally  developed  for  analyzing 
sea  clutter,  its  applicability  to  ground  surface  returns  is  not  fully  understood.  Nevertheless,  it  has 
been  successfully  used  for  describing  a  wide  range  of  land  clutter  {60,  61). 

The  distributions  of  the  vv-polarized  ground  response  for  the  previously  imaged  9  m  x  19  m 
scene  are  displayed  in  figure  14.  Four  ground  roughness  levels  are  considered:  hrms  =1.2  cm, 

1.6  cm,  2.0  cm,  and  2.4  cm;  lc  is  kept  constant  at  14.93  cm  for  the  four  cases.  (The  images  for 
hrms  =  1.2  cm  and  1.6  cm  are  shown  in  figures  11a  and  1  lb,  respectively.)  The  histograms  are 
formed  directly  using  the  image  data  for  the  scene  simulated  with  the  targets  removed.  The 
distribution  p(crcdB  \  is  seen  to  provide  an  excellent  fit  to  the  full-wave  simulated  response: 

(a)  v=  100,  £=3.31  x  10“5;  (b)  v  =  12,  £  =  1.25  x  10“4;  (c)  v  =  10,  £  =  1.65  x  10“4;  and 
(d)  v  =  8,  £  =  2.06  x  10  4.  A  smaller  value  for  v  signifies  a  longer  -tail”  in  the  amplitude 
distribution  of  the  ground  response — implying  there  is  a  greater  probability  for  the  occurrence  of 
strong  scattering  points  within  the  scene.  A  manifestation  of  this  effect  is  readily  seen  in  the 
images.  For  example,  even  though  the  mean  value  of  the  ground  response  for  hrms  =  2.4  cm  (case 
(d))  is  only  slightly  greater  than  that  for  hrms  =  2.0  cm  (case  (c)),  the  image  for  case  (d)  appears 
more  cluttered  due  to  the  presence  of  a  larger  population  of  bright  clusters  of  scattering  centers. 

Figure  14  shows  that  for  a  relatively  smooth  ground  surface  (case  (a)),  Rayleigh  scattering  is  a 
reasonable  descriptor,  but  as  roughness  increases,  the  scattering  statistics  tend  to  deviate  from 
Rayleigh.  One  of  the  generally  held  tenets  of  the  Rayleigh  model  is  that  the  number  of  scatterers 
within  a  resolution  cell  must  be  large.  As  such,  the  central  limit  theorem  can  be  invoked  in 
characterizing  the  received  signal  as  a  vector  sum  of  many  samples  of  random  phasors.  Given 
that  requirement,  as  well  as  the  high  resolution  of  the  images  considered  here  (i.e.,  few  scatterers 
in  each  resolution  cell),  the  observation  that  the  Rayleigh  model  can  be  applied  for  small 
roughness  scales  is  rather  unexpected.  However,  it  is  conjectured  here  that  because  of  the  wide 
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observation  angle  (in  </>)  and  the  multi-aperture  (in  6)  configuration  used  for  the  imaging 
geometry,  the  received  signal  for  a  gently  varying  surface  can  still  be  thought  of  as  the  sum  of 
many  random  independent  samples  of  the  scattered  field. 

The  deviation  from  the  Rayleigh  scattering  process  observed  for  cases  (b)-(d)  can  be  attributed 
to  the  appearance  of  an  increased  number  of  strong  point-like  scatterers  in  a  rougher  surface. 
These  point-like  scatterers,  or  scattering  centers,  tend  to  radiate  in  a  more  directional  and  less- 
random  manner.  Thus,  the  requirement  for  independent  random  sampling  may  be  violated.  In 
addition,  note  that  surface  self-shadowing  and  higher  order  (or  multiple)  scattering  begin  to  play 
a  more  significant  role  in  determining  the  total  signal  return  at  larger  surface  roughnesses,  and 
these  effects  are  further  accentuated  at  low-grazing  angles.  That  confirms  the  well-known  fact 
that  shadowing  and  multipath  effects  can  lead  to  a  departure  from  Rayleigh  scattering  statistics 
(62).  Similar  analysis  as  above  demonstrates  that  expression  13  is  also  a  reasonable  model  for  the 
ground  response  in  the  hh-,  hv-,  and  v’/i-polarized  images.  Of  the  four  polarization  combinations, 
the  distributions  for  hh  are  seen  to  deviate  from  Rayleigh  statistics  the  most,  while  the 
distributions  for  the  cross-polarized  components  are  seen  to  lie  somewhere  between  those  of  the 
vv  and  hh  cases.  This  observation  is  consistent  with  previously  reported  results  on  the 
dependency  of  scattering  distribution  on  polarization  (63). 
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Figure  14.  Probability  density  of  vv-polarized  ground  response:  histogram  shows  the  distribution 
inferred  from  simulated  images;  solid  line  is  distribution  calculated  using  expression  13. 
Ground  with  randomly  rough  surface:  (a)  hrms  =1.2  cm,  lc  =  14.93  cm;  (b)  hrms  =1.6  cm, 
lc  =  14.93  cm;  (c)  hrms  =  2.0  cm,  lc  =  14.93  cm;  and  (d)  hrms  =  2.4  cm,  lc  =  14.93  cm.  Ground 
electrical  properties:  sr  =  8,  ad=  10  mS/m. 
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5.  Conclusion 


A  FDTD-based  solver  is  proposed  for  characterizing  the  electromagnetic  scattering  behavior  of 
rough  ground  surfaces  at  low  depression  angles.  The  simulation  approach  employs  a  finite-extent 
surface  perturbation  model  for  the  interface  and  a  coherent  error  subtraction  technique  for 
reducing  numerical  dispersion  artifacts.  The  regions  of  validity  of  the  solver  have  been 
investigated  for  1-D  surfaces  by  comparing  Monte  Carlo  scattering  results  to  those  from  a  high- 
order  SIE  approach  for  various  surface  parameters  and  incidence  angles.  Backscattering 
coefficients  of  2-D  random  surfaces  as  computed  with  the  MPI-based,  parallelized  FDTD  code 
display  good  agreement  with  measurement  data  available  in  literature.  The  study  indicates  that 
the  overall  FDTD  framework  described  herein  is  a  promising  approach  for  modeling  the  grazing 
scattering  effects  of  a  wide  range  of  surfaces  encountered  in  practical  radar  sensing.  Inaccuracies 
have  been  noted,  however,  when  the  surfaces  exhibited  at  once  small  rms  slope  and  large 
correlation  length — parameters  yielding  geometrical  variations  that  are  difficult  to  model 
correctly  within  the  stair-stepped  approximation  of  the  FDTD  grid.  Although  not  shown  here,  if 
the  results  were  extrapolated  for  the  treatment  of  surfaces  exhibiting  large  rms  slope  and  small 
correlation  length,  the  FDTD  algorithm  is  expected  to  encounter  difficulties  for  that  case  as  well. 
The  aforementioned  limitations  and  errors  can  be  mitigated  by  refining  the  computational  lattice. 
Another  possible  improvement  would  be  implementing  a  non-uniform  lattice  that  uses  finer  cells 
around  the  ground  interface  region  (64),  or  using  stretched  Yee  cells,  in  which  one  dimension  is 
larger  than  the  others  (35).  With  respect  to  the  proper  selection  of  the  rough  surface  spectra  (or 
autocorrelation  function),  this  work  focuses  mostly  on  Gaussian  surfaces.  However,  the 
numerical  analysis  methods  are  more  general  and  can  be  applied  to  realistic  surfaces  of  arbitrary 
spectra  (including  multi-scale  models). 

As  compared  to  conventional  SIE,  the  proposed  FDTD  algorithm  is  expected  to  be  a  more 
computationally  viable  option  for  analyzing  grazing-angle  scattering  phenomena,  since  it 
necessitates  a  smaller  simulation  domain  with  horizontal  dimensions  that  are  independent  of  the 
incidence  angle.  Additionally,  the  FDTD  technique  can  easily  accommodate  hybrid  scattering 
scenarios,  where  targets  of  almost  arbitrary  shape  and  material  composition  can  be  included 
together  with  a  rough  air-ground  interface,  as  well  as  soil  inhomogeneities. 

In  demonstrating  the  usefulness  of  the  developed  solver  as  pertinent  to  forward-looking  radar 
sensing,  the  effects  of  surface  clutter  on  multistatic  target  imaging  are  illustrated  with  large-scale 
simulations  of  a  realistic  scene  consisting  of  targets  embedded  in  a  rough  ground  environment.  A 
study  of  the  ground  response  in  the  image  domain  shows  that  the  ^-distribution  provides  a 
reasonable  model  for  characterizing  the  statistics  of  the  rough  surface  multistatic,  wideband  radar 
return. 
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Collectively,  the  formalisms,  results,  and  observations  featured  in  this  study  represent  the  first 
step  taken  in  understanding  and  predicting  the  impact  of  ground  surface  clutter  on  the  detection 
performance  of  a  low-frequency  imaging  radar  operating  at  close-to-grazing  angles.  Future  work 
will  include  more  realistic  models  of  the  SIRE  radar  geometry  (which  could  better  be  described 
by  a  near-field  configuration)  and  its  antenna  patterns,  as  well  as  the  roughness  characterizing 
the  terrain. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


APSWF 

approximate  prolate  spheroidal  wave  function 

ARL 

U.S.  Army  Research  Laboratory 

CPU 

central  processing  unit 

DSRC 

DoD  Supercomputing  Resource  Center 

FDTD 

finite-difference  time-domain 

FMM 

fast  multipole  method 

hh 

horizontal-horizontal 

hv 

horizontal-vertical 

KA 

Kirchhoff  approximation 

LCN 

locally  corrected  Nystrom 

MoM 

method  of  moments 

MPI 

Message-Passing  Interface 

PMCHWT 

Poggio,  Miller,  Chang,  Harrington,  Wu,  and  Tsai 

PML 

perfectly  matched  layer 

RCS 

radar  cross  section 

rms 

root  mean  square 

SDFMM 

steepest  descent-fast  multipole  method 

SIE 

surface  integral  equation 

SIRE 

synchronous  impulse  reconstruction 

SPM 

small  perturbation  method 

TF/SF 

total  field/scattered  field 

UHF 

ultra  high  frequency 

uxo 

unexploded  ordnances 

vh 

vertical-horizontal 
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VHF 

vv 


very  high  frequency 
vertical-vertical 
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